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Abstract 



We study fluidized granular gases in a stationary state determined by 

the balance between an external driving and the bulk dissipation. The two 

considered situations are inspired by recent experiments, where the gravity 

C^ . plays a major role as a driving mechanism: in the first case gravity acts 

C . only in one direction and the bottom wall is vibrated, in the second case 

i-^ \ gravity acts in both directions and no vibrating walls are present. Simulations 

C ' performed under the molecular chaos assumption show averaged profiles of 

fi ' density, velocity and granular temperature which are in good agreement with 

the experiments. Moreover we measure the velocity distributions which show 

CN . strong non-Gaussian behavior, as experiments pointed out, but also density 

J^ \ correlations accounting for clustering, at odds with the experimental results. 

f^ ' The hydrodynamics of the first model is discussed and an exact solution is 

^~~^ • found for the density and granular temperature as functions of the distance 

^^ ■ from the vibrating wall. The limitations of such a solution, in particular in a 

O , broad layer near the wall injecting energy, are discussed. 
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I. INTRODUCTION 

In general granular materials [|I|, since the presence of dissipative forces, are not equilib- 
rium systems neither from the configurational point of view nor from the dynamical point 
of view. A statistically stationary state can be produced by the competition between the 
dissipation due to the inelastic collisions among the particles and the energy injection due 
to an external source which prevents the system from cooling and come to rest. 

Usually granular gases are considered in the homogeneous cooling regime, less frequently 
they are studied in a stationary regime where energy flows into the system from some external 
source (stochastic driving, vibrating plates, shear, etc.) and dissipates by means of inelastic 
collisions. A sufficient condition to prevent strong density instabilities (such as those found 
by Du et al. [Q), seems to be the presence of an even minimal, but spread, temperature 
source |^. 

Many evidences, by mean of computer simulations, have been found that different kinds 
of density instabilities, like clustering (density gradients growing on time scales faster than 
typical hydrodynamics scales) or inelastic collapse (the local divergence of collision rate so 
that an infinite number of collisions occurs in a finite time) may emerge in a cooling granular 
assembly, that is a granular gas losing his starting kinetic energy because of dissipative 
collisions. It has also been shown that the velocity distribution of particles in the free cooling 
state with homogeneous density has overpopulated high energy tails ~ exp{—Av) ^j^. 

When granular gases are driven in some way to balance the loss of energy due to collisions, 
a stationary state may be observed. The first model of randomly driven granular gas was 
proposed in P|. It showed pathologies in the density and granular temperature profiles 
but also a breakdown of thermodynamic limit. Another randomly driven model was then 
proposed to offer a different insight into the kinetics of granular gases [^. In this model 
the driving mechanism is a stochastic energy source acting on every particle as a heat 
bath with a fixed temperature Tp and a fixed viscous damping with characteristic time 
T. In the stationary "collisional" regime (characterized by a collision time much lower 
than r) the gas showed a fractal distribution of density and a distribution of velocities 
with overpopulated (non- Gaussian) high energy tails. The homogeneous solution of the 
corresponding Boltzmann-Enskog equation has been analytically studied ||^ showing that 
~ exp{—Av^''^) high energy tails are expected. 

The aim of this work is to study a class of models for driven granular gases where the 
efficiency of the energy injection is guaranteed by the presence of gravity, taking inspiration 
by some recent experiments @,0: in these experiments a bottom confining wall is the 
source of granular temperature while gravity forces the particles to return in contact with 
this source. We are interested in very diluted systems, where the granular material behaves 
as an inelastic gas, rather than dense granular fiows, where many static effects, such as 
clogging, arching or bubbling appear. Such systems have been studied in relation with 
compaction dynamics or slow dense chute flows 0. The study is based on Direct Simulation 
Monte Carlo, but we also discuss (for one of the models) the hydrodynamic theory. The 
flrst version of the model (gravity in only one direction and vibrating bottom wall) has 
been previously studied in the one- dimensional case, that is a vibrated column of grains 
under the force of gravity [|1T| and the transition or the coexistence of different phases (gas, 
partially fluidized and condensed) was investigated. In two dimensions experiments [|T2|, 



simulations |T3| and theories |]14[ have analyzed a vertical system of grains with gravity 
and a vibrating bottom wall (with different kinds of vibration) searching for a simple scaling 
relation between global variables as the global granular temperature Tq or the center of mass 
height hem as function of the size of the system A^, the typical velocity of the vibrating wall 
V or the restitution coefficient r. In all these calculations the authors did not pay too much 
attention to the hydrodynamic profiles of the system, always assuming a constant granular 
temperature ( "isotherm atmosphere" ) and a density profile exponentially decaying with the 
height, as in the case of a Boltzmann elastic gas under gravity. One of the results of this 
work, discussed in section V, is that also in the dilute regime that one can study by means 
of Monte Carlo methods the use of these assumptions is not obvious, in particular when 
trying to solve the global balance between external energy injection and bulk dissipation 
due to inelastic collisions among particles. It must be also noted that the general validity of 
a hydrodynamical description is still subject of debate in the case of granular gases far from 



the elastic limit (a review of hydrodynamical problems is |T3]). 

In section II we present the two versions of the model, in section III and IV we illustrate 
the results, in section V we discuss the hydrodynamics of the model in its first version, and 
finally we draw the conclusions. For the sake of completeness and in order to make the paper 
self-contained we included in Appendix A a brief description of the Direct Simulation Monte 
Carlo of the Boltzmann equation and in Appendix B the expressions of the dimensionless 
coefficients appearing in the hydrodynamic equations of section V. 

II. THE MODELS 

We introduce two bi-dimensional models both consisting of A^ identical smooth disks 
of diameter a and mass m = 1 subject to binary instantaneous inelastic collisions which 
conserve the total momentum 

v'l + V2 = Vi + V2 (1) 

and reduce the normal component of the relative velocity 

K - v^) ■ fi = -r((vi - V2) ■ n) (2) 

where r is the normal restitution coefficient (r = 1 in the completely elastic case) and 
n = (xi — X2)/cr is the unit vector along the line of centers xi and X2 of the colliding disks 
at contact. With these rules satisfied the post-collisional velocities are: 

1 + r 
vi' = vi ^((^1 " ^2) ■ n)n 

1 + r 
V2 = V2 + ^—((vi -V2) -11)11 (3) 

In addition, the particles experience the external gravitational field and the presence 
of confining walls. With respect to previous works [Q the energy necessary to prevent the 
cooling of the system due to the inelastic collisions is not provided by a heat bath: in the 
present paper the energy feeding mechanism is of two types according to the two numerical 
experiments we perform. 



In model a), illustrated in fig. (0) and inspired to a recent laboratory experiment [0 and 
a numerical experiment [^ , the "apparatus" consists of a plane of dimension L^, x Ly 
inclined by an angle 9 with respect to the horizontal. The particles are constrained to 
move in such a plane under the action of an effective gravitational force Qe = g sin 6 
pointing downward. In the horizontal direction there are periodic boundary conditions. 
Vertically the particles are confined by walls, both inelastic with a restitution coefficient 
Tw (difference of restitution coefficients for particle-particle interactions and particle- 
wall interactions will be discussed below). Besides, the bottom wall vibrates and 
therefore injects energy and momentum into the system. The vibration can have 
either a periodic character (as in 0) or a stochastic behavior with thermal properties 
|16|). In the periodic case, the wall oscillates vertically with the law F^ = 



as m 



Au, sm{ujujt) and the particles collide with it as with a body of infinite mass, so that 
the vertical component of their velocity after the collision is v'y = —r^Vy + (1 -l- r^)K; 
where V^ = A^ui^ cos(a;^t) the velocity of the vibrating wall. In the stochastic case we 
assume that the vibration amplitude is negligible and that the particle colliding with 
the wall have, after the collision, new random velocity components Vx G {~oo, +cxd) 
and Vy G (0, +cx3) with the following probability distributions: 



P{vy) = l^exp(--i) (4) 



In model b) sketched in fig. (0) the "set-up" is a two dimensional channel of depth Ly 
and of length L^., vertically confined by a bottom and a top inelastic wall, with periodic 
boundary conditions in the direction parallel to the flow. The channel is tilted up by an 
angle with respect to the horizontal so that gravity has both components Qx = 9 sin (p 



and Qy = g cos 0. This model mimics the experiment performed by Azanza et al. []10 
where a stationary flow in a two dimensional inclined channel was observed at a point 
far from the source of the granular material. The assumption of periodic boundary 
conditions in the direction of the flow is consistent with the observed stationarity, due 
to the balance between the gravity drift and the damping effect of inelastic collisions 
(for a discussion of the possible regimes that can be shown by one particle in presence 
of this balance, see ii^ 



The chosen collision rule excludes the presence of tangential forces, and hence the rota- 
tional degrees of freedom do not contribute to the description of the dynamics. 

Under the assumption of molecular chaos, stating that P2(x, x -|- an, vi,V2,t) = 
P(x, vi,t)P(x + an, V2,t) where P2 and P are the probability density functions for two 
particles and one particle respectively, it is possible to write down the Boltzmann equation 
(eq. (y) in section V), which can be solved by means of Monte Carlo methods. Here we 
used a simplified (but still efficient) version of the Direct Simulation Monte Carlo scheme 
proposed by Bird [O]. With respect to the original version of the algorithm the clock which 



determines the collision rate is replaced by an a-priori fixed collision rate via a constant 
collision probability Pc given to every disk at every time-step At of the simulation, in such 
a way that the single-particle collision rate is x ~ Pc/^t. The colliding particle then seeks 
its collision partner among the other particles in a neighborhood of radius r^, choosing it 
randomly with a probability proportional to their relative velocities. Moreover in this ap- 
proximation the diameter a is no more explicitly relevant but it is directly related to the 
choices of pc and rs in a non trivial way: in fact the Bird algorithm allows the particles to 
pass through each others, so that a precise diameter cannot be defined and estimated as a 
function of Pc and tb- The Bird scheme is described in more details in the Appendix A. 

The agreement between our simulations and the inspiring experiments, justifies the sim- 
plifying assumptions considered for our model, i.e. assuming molecular chaos and neglecting 
tangential forces. Nevertheless, as a partial check, we try a modified version of model b) 
where the tangential forces may affect the post-collisional velocities of the particles. As re- 
ported below, the introduction of such forces does not change the behavior of the measured 
quantities. 

III. DISCUSSION OF THE RESULTS: 
MODEL A) 

Simulations of the first model, the inclined plane with a bottom wall injecting energy, 
have been performed for different choices of the number of disks N, the normal restitution 
coefficient r, the dimensionless width of the plane N^j = L^/vb and the parameter measuring 
the rate of energy injection from the wall, that is the temperature T^ in the stochastic case 
and the amplitude and frequency A^;, u^ in the periodic case. 

Let us show how numerical simulations with the molecular chaos assumption reproduce 
the main results obtained in experiments [§,0] and in high performance computer simula- 
tions [jl6| of inelastic hard disks. 



Snapshots of the systems and time-averaged density profiles are shown in fig. (|) for the 
case of randomly vibrated wall. We are in the presence of a highly fluidized phase of the 



type Isobe and Nakanishi []T6[ call granular turbulent: looking at the time evolution of the 
density distribution of the system and of the coarse-grained velocity field one observes an 
intermittency-like behavior with rapid and strong fiuctuations of the density, sudden explo- 
sions followed by large clusters of particles traveling downward, coherently, under the action 
of gravity. Of course more dense and ordered phases (that one can expect at lower values 
of energy injection) are not reproducible with the Direct Simulation Monte Carlo, as strong 
excluded volume effects appear and the assumption of negligible short range correlations 
fails. 

In the figures (^ and (|5|) we display the horizontal velocity distributions, for the stochas- 
tic case. In fig. (^) distributions for different T^ are shown: the data collapse is obtained 
by rescaling the velocities by ^/T^. Instead, in fig. (H) we show the velocity distributions 
of particles contained in stripes at different heights from the wall, again rescaled by y^T(jj) 
(their own variance) in order to obtain the data collapse. It appears that the distributions 
are non-Gaussian and their broadening (that is the granular temperature T{y)) is density 
dependent. This dependence is shown in fig. (^) as well as its dependence upon the height. 
An analogous dependence has been shown in references [§] where the granular gas was driven 



by a homogeneous heat bath, showing a power law T ~ n~^ with (5 ~ 0.8, while in this case 
it seems (5 ~ 0.88. 

The case of periodically vibrated wall is illustrated in figures (^ and (^. One can see the 
density profiles (together with a snapshot of the system) and the distribution of horizontal 
velocities in two different regimes: for g^ = —1 a. non-Gaussian distribution is obtained, 
while a distribution close to a Gaussian appears when g^ = —100. This trend towards 
a Gaussian, as the angle of inclination is raised up, reproduces exactly the experimental 
observation of KudroUi and Henry |^ (where the angle of inclination of the plane was raised 
up from 6 = 0.1° to ^ = 10°) and can be explained as an effect of the increase of the collision 
rate with the wall which "randomizes" the velocities in a more efficient way: this resembles 
the heath bath model [^ where one passes from the non-Gaussian regime to the Gaussian 
one increasing the ratio between the heating rate and the collision rate. 

In order to characterize the spatial clustering we have studied the cumulated particle- 
particle correlation function: 



CBiyAy)it,R) = ^ r^. TT Yl e(i?-|x,(t)-x,(t)|) (6) 

where B{y, Ay) is a horizontal stripe contained between y + Ay/2 and y — Ay/2. After 
having checked that the system has reached a stationary regime, we have computed the 
time-average of the correlation function, that is 

1 [^ 

J- -to J to 

which is independent on time if T » to- In the figure (P) we show the C{R) vs. R for 
different stripes B{y, Ay). We observe a power law behavior 

CBiyAy)iR) - R''^'^ (8) 

In the case of homogeneous density d2 is expected to be the topological dimension of the 
stripe, that is d2 = I ii R^ Ay and ^2 = 2 if i? <^ Ay. 

Clustering, whose signature is a value of the correlation dimension d2 lower than the 
topological dimension, appears in the stripes with not too high densities, where an exponent 
smaller than 1 is measured (the fit is performed in th region R ^ Ay). The evidence of 
clustering is at odds with the observation of Kudrolli and Henry 0: they report, in fact, the 
absence of clustering by measuring the distribution of the number of particles in boxes of 
fixed dimensions spread all over the inclined plane. This observation is perhaps due to the 
fact that in the statistical analysis employed in ref. the number of particles in each box is 
considered disregarding their heights, that is they may belong to regions of different densities: 
in such a way the slow decaying tails, expected for the clusterized distributions of the stripes 
at lower densities, are partially hidden by the Poissonian (homogeneous) distribution of the 
stripes at higher density. Moreover, even from the global density distribution measured in 
their work, a tail decaying slower than a Poissonian cannot be clearly ruled out. 
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IV. DISCUSSION OF THE RESULTS: 
MODEL B) 

Let us now show the results for the second model, the inclined bidimensional channel. 

In figures (^) and (pT]) the hydrodynamical fields n{y) (number density), Vx{y) (velocity 
component parallel to the fiow), T{y) (granular temperature) are shown as functions of 
the distance from the bottom wall y. The velocity, the temperature and the height are 
made dimensionless by rescaling them by y/gxfB, Qx'^b and r^ respectively. The profiles 
well reproduce those measured experimentally by Azanza et al. [T^: they show a critical 



height H of about six times the radius tb which corresponds to the separation between two 
different regimes of the cooling rate. In a mean field framework the local rate of dissipation 
due to the inelastic collisions (as already stated before) is C oc nT^^"^. 

This can be understood simply noting that the collision rate is proportional to the local 
density and to the local relative velocity of the particles (v^), while the change in the 
granular temperature induced by every collision is proportional to the temperature T. The 
quantity C, = riT^/'^ as a function of y is shown in fig. (^). The cooling rate decreases 
exponentially and is reduced under 1/100 of its maximum value at about the observed 
critical height H ^ Qtb, accounting for the difference between a collisional regime and a 
ballistic one. 

With respect to the velocity and temperature profiles in figure (|T0|), we note here that 
quite unphysical features appear: in particular the quite strong slipping effect near the 
bottom wall is in contrast with the experimental findings. We think that this is due to a 
wrong modelling of the particle-wall collision events. 

The restitution coefficient used in our model has to be considered as an effective param- 
eter describing the energetics of collisions. It should depend on the details of the colhsion 
event, in principle even on the relative velocities of the colliding particles. In the experiment 
the bottom wall was covered with particles identical to the flowing ones with a spacing 
bounded between and 0.8 mm: however the particles are stuck to the bottom wall so that 
the collision event is completely different from a two-particles collision. 

Using a lower effective restitution coefficient for the wall r^ (see figure (PUJ)) we obtain 
a better agreement with the experimental profiles. In particular, both temperature and 
velocity profiles seems to go to zero near the bottom, although we cannot really rule out 
slipping effects {vx{y = 0) 7^ 0). 

We have also studied the distribution of horizontal velocities in stripes at different heights 
(here the mean values are height dependent). These are displayed in fig. (|1^) showing the 
emergence of a non-Gaussian behavior mainly in the case with r^ < r and only in the stripes 



near the bottom wall. The authors of the experiment of ref. |T^ claim that the distributions 
of velocity are very close to the Gaussian and try to fit their data with the rheological model 
proposed by Jenkins and Richman |T^, which postulate a quasi-Gaussian equilibrium to 



calculate the transport coefficients. Near the bottom wall Gaussian approximation is far 
from obvious, as shown by the results of our simulations: this is an effect of the inelasticity 
of the collisions but also of the proximity of the boundary. 

Finally we have investigated the homogeneity of the density: the figure ([I^) shows the 



previously defined function CB{y,Ay)iR) for stripes at different density. It appears again a 
clustering effect, with a correlation dimension ranging from 1 (homogeneous stripes) to 0.2 



(highly clusterized stripes). In the figure it is also shown the very small distance region, 
R < rB, where homogeneity should be recovered. Since in our simulation, Ay ^ tb, we 
expect d{y) = 2 in this region. 

We consider the comparison between our simplified model and the experimental profiles 
quite satisfactory: this seems to suggest that introducing further physical details should be 
irrelevant at this description level. However we briefiy report the results obtained with a 
slightly modified version of the model, including the effects of tangential forces. Such forces 
play a key role in dense granular flows P,pO|, being responsible for arching. On the other 



hand the present results suggest that in the case of diluted systems they act similarly to the 
normal forces without introducing noticeable effects. 

The introduction of tangential forces in the model studied accounts for a new collision 
rule: 

(v'l - V2) • n = -r"((vi - V2) ■ n) 
(v'l - V2) ■ t = -r*((vi - V2) • t) 

where we replace the single restitution coefficient with a pair of parameters r" and r*, 
respectively due to the effect of normal and tangential collision forces {t is a unit vector 
perpendicular to n) . Analogously, the restitution coefficient r^ splits in two new parameters 
r^ and r^. The results of simulations with several choices of the enlarged set of parameters 
do not show qualitative differences: setting tangential restitution coefficients lower than one 
is equivalent to enhance the dissipation in the original model. 

Just as an example of this, we show figure (|lD, where the extremal case of a vanishing 
tangential restitution coefficient is reported. Note that the profiles are similar to those shown 
in figure ( |ri|) where a low r^j = 0.4 was used. 

V. DISCUSSION OF THE HYDRODYNAMICS: RESULTS AND PROBLEMS 

The Boltzmann equation for the two models introduced in this paper (in two dimensions) 
reads: 

(I +V.V + .. A) ;(.,.,,).,(;,;) (9) 

Jif, f) = ajdY,j dhe{h ■ v,)(n ■ V,,) [r-V(x, v', t)/(x, v'^, t) - /(x, v, t)/(x, Vi, t)] 

(10) 

Here n is the unit vector along the line joining the centers of the colliding particles at 
contact, Vr = V — vi is the relative velocity of the colliding disks, is the Heaviside step 
function, v' and v'^^ are the precoUisional velocities leading after collision to velocities v,vi. 

The equation (||) must be completed with the boundary conditions in order to describe 
the microscopic evolution of the whole system. 

The difficulty of solving the Boltzmann equation (H) can be bypassed substituting the 
microscopic description given by /(x, v,t) with the averaged macroscopic description given 
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by the hydrodynamic fields: the number density field n{'x,t), the velocity field v(x, t) and 
the granular temperature field T(x, t). These quantities are given by 



n X, 



t)= j dwf{^,^,t) (11) 



u(x, t) = -^ / civv/(x, V, t) (12) 



kBT{^, t) = — — / rfv^ ^^^/ X, V, t) 13 

n{x,t) J 2 

Multiplying the Boltzmann equation (O) by 1 or v or m(v — u(x, t))^/2 and integrating 



over vi one can derive l2Jj2^ the equations of fiuid dynamics: 



Dr) 

— + nd,Ui = (14) 



mn—^ = -djTij + ngiTTi {i = 1, 2, 3) (15) 



n = -diQi - TijdjUi - CnksT (16) 

where di = d/dxi (for the sake of compactness we use here the notation x — *> xi and y — >■ 
X2) and D/Dt = d/dt + u ■ V is the Lagrangian derivative, e.g.: -^F^x, t) = J^F(0(xo, t),t) 
with (j)['XQ,t) the evolution after a time t of xq under the velocity field u. In the above 
equations 

r.. = /*,„(.. -,..)(..-«.)/(x.v,«) (17) 

is the stress tensor, g is the volume external force (gravity in our case), 

/in 
C?Vyfi|v-u|V(x,V,t) (18) 

is the heat fiux vector and 

^^""'^^^ 8r(5/2)nfcBr y ^^ly ^^2|vi-V2|V(x,Vi,t)/(x,V2,t) (19) 



is the cooling rate due to dissipative collisions. 

The set of equations (|14|) -(|T^) become closed hydrodynamic equations for the fields n, u 
and T when Pij, q and ( are expressed as functionals of these fields. This is obtained, for 
example, expressing the space and time dependence of / in terms of the hydrodynamic fields 
and then expanding / to first order (the so-called Navier-Stokes order) in their gradients, 
with the exception of ( which requires an expression of / to the second order of gradients to 
be consistent with the other terms. With this approximation the equations (|T^)-(^) include 
the contributions up to the second order in the gradients of the fields. 

Calculations of the closure of the hydrodynamic equations for granular media have been 
performed with some approximations restricting the validity of the results to the low dis- 
sipation or quasi-elastic limit [|^. More recently |2l[] the analysis has been extended to 



arbitrary inelasticity giving closed expressions for the momentum and heat fluxes and for 
the cooling rate (. 

We follow these more recent results |]2T]] and write down the hydrodynamics for the 
model a) presented in this paper (gravity in one direction and vibrating bottom wall, i.e. 
g = {0,ge) and Qe < ), with the following assumptions: the fields do not depend upon x 
(the coordinate parallel to the bottom wall), i.e. d/dx = 0, and the system is in a steady 
state, i.e. d/dt = 0. The continuity equation (|T^) then reads ^{n{y)uy{y)) = and this can 
be compatible with the bottom and top walls (where nVy = 0) only if n{y)vy{y) = 0, that is 
in the absence of macroscopic vertical flow. The equations are written for the dimensionless 
fields T = ksT /{—gema) and ft = na'^, while the position y is made dimensionless using 
y = y/a. Finally for the pressure we put p{y) = T22 = n{y)kBT{y). With the assumption 
discussed above the equations of Brey et al. pl[ read: 



^(niy)f{y)) = -n(y) (20) 



^ '^.Qr{y) + C{r)n{y)f{yf' = (21) 



h{y) dy 
where Qr{y) is the granular heat flux expressed by 

Qr{y) = A{r)f{yy/'-f{y) + B{r)^^-n{y) (22) 

In the above equations A{r), B{r) and C{r) are dimensionless monotone coefficients, 
all with the same sign, explicitly given in the Appendix B. In particular -B(l) = and 
C(l) = 0, i.e. in the elastic limit there is no dissipation and the heat transport is due only 
to the temperature gradients, while when r < 1 a term dependent upon ^ ln(n(^)) appears 
in Qr{y)- The use of dimensionless fields eliminates the explicit g dependence from the 
equations, that remains hidden in their structure (the right hand term of equation ^, that 
is due to the gravitational pressure gradient, disappears in the equation for g = Gi). 

A change of coordinate can be applied to eqs. (^),(^) in order to obtain a simpler 
form: 
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y ^ l{^y) = / n^y')dy' (23) 





It follows that when y spans the range [0, Ly], the coordinate / spans the range [0, cr/Lx]. 
With this change of coordinate it happens that 

and the first equation ( pO] ) reads: 

|(n(/)t(/)) = -1 (25) 

from which is immediate to see that 



H = n{l)f{l) + / (26) 

is a constant, i.e. ^^H = 0. This is equivalent to observe that 

Piy)-9 rniy')dy' (27) 

Jo 

is constant which is nothing more than the Bernoulli theorem for a fluid in the gravita- 
tional field with the density depending upon the height. 



The relation (p6|) is verified by the model simulated in this work in the figure (|1^) for 
almost all the height of the container, apart of the boundary layer near the bottom driving 
wall. 

Using the coordinate I introduced in (^) and the elimination of n{l) using the recognized 
constant, that is 

n{l) = ^Lzl (28) 

the second equation (|2l|), after some simplifications, and after a second change of coor- 
dinate / -^ s{l) = H — I, becomes: 

a(r)s (P ~, , a(r)s f d ~, .\ Bir) d ~, . ~, ,i/o , , 



f{sy/^ds^ '' 2f{sy/^ yds '7 f(s)i/2ds 

where a(r) = {A{r) — B{r))/C{r-) =, (3{r-) = {A{r) — ^B{r))/{C{r)) are numerically 
checked to be positive (see Appendix B) for values of r not too low (about r > 0.3) and are 
divergent in the limit r — ;► 1. 
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The equation (^) become a linear equation in T{s) as soon as the change of variable 
z{s) = T{sY/'^ is performed: 



(f d 

2a{r)s—^z{s) - 2/?(r)— ^(s) + szis) = (30) 

as"^ as 



giving the solution: 



z{s) = As"' J,,(,)(/3'(r)s) + 5s"'MiV„,(,)(/5'(r)s) (31) 

where Ja' and Na' are the Bessel functions of the first kind and the second kind respec- 
tively, a'{r) = (a(r) +/3(r))/(2a(r)) is real and positive, (3'{r) = (l/(2a(r)))"^/^ is real and is 
considered in its positive determination, moreover they present the elastic values a'(l) = 1 
and f3'{r — ;> 1) = (see appendix B), while A and B are constants that must be determined 
with assigning the boundary conditions. 

Then we can derive the expressions for T{1) and n{l): 

f{l) = {H- lf'^'^^\AJ^,^,^{(3\r){H - 0) + 5iV,,(,)(/3'(r)(iJ - l))f (32) 



""^ ^ " {AJ^'(r){P'{T){H - I)) + BN^,^r){P'{T){H - l))Y ^'^^> 

To calculate the expressions of T and n as a function of the original coordinate y one 
needs to solve the equation 

putting in it the solution (|33D . However one can obtain a comparison with the numerical 
simulations using the new coordinate /. The main problem, at this point, is a discussion of 
the boundary conditions needed to eliminate the constants if, A and B. 

One could impose that n{lmax) = at Imax = o"/-^x- From this condition immediately 
follows that H = cx/Lx. A second condition can be obtained imposing a vanishing derivative 
of the temperature at Imax, that is 

The third condition is the most delicate: it must contain the rate of energy injection 
coming from the vibrating wall. This rate depends upon the parameter T^ (or A^ and 
uj^) and upon the particles flux impinging on the wall $ = n^L^/U^ where n+ is the number 
density of particles approaching the wall and W^ is their velocity averaged near the wall. The 
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first may be simply estimated as n+ = n(0)/2. Moreover, if the velocity of the macroscopic 
fiow is zero, the average velocity of the impinging particles is due only to fiuctuations of u, 
that is I)^ ~ A/fc_BT(0)/m. In a collision with the wall, the average energy gain is given by: 



— v'r — vM = — 
2 VI I I I ; 2 



AE^ = -(|v'|2 - |vP) = -[«)2 + {v'^y - {v^y - [vyy] (36) 



which is different for the stochastic or the periodic case, respectively: 

a:^;;: = ^^ - ksno) (37) 



AE^, = ^ ' 2 " " - ^ 2 ^^^^ 

obtained straightforwardly from the eq. ( ^6l) assuming no correlations between the ve- 
locity of the wall and that of the approaching particles. 

Then a non-closed expression for the rate of energy injection coming from the wall reads: 



3, ^ , ,,^ kBT{0) L,n(0)(A;BT(0))3/2 



W^, = AE^,<^ = -ksT^LMO)\/^^- ' ' 'l^' " (39) 

m(l + r^)2A>% ,^, [k^m (1-^'). ,,,(A:bT(0))3/2 



W^p = AE^p$ = ^ ' ;^ " L^n{S))\^^^^ - ' ^ ""' LMO y ''' (40) 



The above expressions are useful to establish the third needed boundary conditions. In 
order to do that, they must be compared with the energy dissipation rate due to inelastic 
collisions. The local dissipation rate is given by ({y)kBT(y) = C {r)an{y){kBT{y)Y^'^ / y/m 
(see Appendix B). The instantaneous balance between energy injection and dissipation in 
collisions reads, then, 

W^-l-r i. r- ,yCt.ny) ^ t^^ip^ r^ ,Mn,rr^ (41) 

J^xJ^y Jo Jo ^y Jo 

where W is W^s or Wu;p. 

Apart of the difficulty of solving the boundary conditions to give an expression of A and 
B as functions of the parameters of the model, one must observe that the hydrodynamic 
description given here is ill-posed from the beginning for what concerns a broad boundary 
layer near the bottom wall. A simple look to the profiles of Uy{y) = Uyijj) / ^J—geTB and 
T{jj) = Ty/^—geTB) (the overlined quantities n, u, T and y are analogue of the dimensionless 
variables n, u, T and y with the assumption rra = 1, fc^ = 1 and a = tb) in the figure (0) can 
give the idea. We expect from the continuity equation (p^, as discussed above, Uyijj) = 
for every y, while a broad region appears with a non-constant and non-monotonic behavior. 
Moreover, even the profile oiT{y) shows an extremal point, in this case a minimum: but from 
equation (^2]) , taking in account the substitution (pSD, it can be seen that the imposition 
■jzT = gives the following relation: 
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where the fraction C /{B — A) = — 1/a is numerically checked to be negative from a value 



of r lower than 0.4 (see appendix B). The relation (42) states that if T > a minimum (that 
is a positive value of -^T) cannot be expected. Similar profiles for T(y), with a minimum, 



have been obtained in other simulations 22 



A tentative fit is presented in figure (18). Here we used three boundary conditions 



obtained directly from the simulations: a value of riT at a certain height yi to obtain 
directly H, the value of T and the value of its derivative at heights 1/2 and y^ respectively, 
with all yi, y2, ys not far from the top wall. In this tentative fit the problems discussed 
above appear clearly: there is a broad region near the bottom wall (see also figure ([T7|)) 
where the theoretical solution of the equations (^)-(^) is qualitatively different from the 
simulation data. 

The qualitative inconsistencies between the observed profiles and the hydrodynamics 
in a broad boundary layer near the vibrating wall are probably due to the high density 
gradients present in this region. The high density gradients represent a numerical but also 
a conceptual problem: it is numerical because the profiles shown in figure (|17D are obtained 
by means of a coarse graining in horizontal stripes B{y, Ay) and so they can be compared 
to the theoretical profiles only if the density in these stripes is approximately homogeneous; 
it is conceptual because this hydrodynamic description is based upon the Navier-Stokes 
approximation, that is an expansion of /(x, v,t) = /[v|n(x, t), u(x, t),T(x, t)] up to the 
first order in the gradients of the fields n, u, T. 

It must be stressed the fact that this boundary layer problem affects the description of 
the whole system in a strong way, as its global behavior (for example the scaling laws for 
the global temperature or the center of mass height, extensively investigated in ||Tl|-p!^) 



emerges from the balance between the bulk dissipation and the injection rate which cannot 
be determined, even qualitatively, by a hydrodynamic study at the level proposed in this 
paper. 

VI. CONCLUSIONS 

We have studied, by means of a Direct Simulation Monte Carlo algorithm, a model of 
granular flow in two different bi-dimensional setups: the first version consists of an inclined 
plane with periodic horizontal boundary condition, a top inelastic wall and a vibrating 
bottom inelastic wall while gravity acts in the direction y perpendicular to the vibrating 
wall and pointing toward it; in the second version gravity acts in both x and y directions 
and the bottom wall doesn't vibrate, therefore resembling a stationary flow along a bi- 
dimensional channel. In both versions of the model we have found a good agreement with 
the analogous experiments [§,0. In particular, the model with the vibrating wall shows 



strong non-Gaussian behavior of the velocities which turns to a Gaussian behavior if the 
angle of inclination is raised up (this should be an effect of the increase of the heating 
rate, as the particles are more frequently in contact with the vibrating wall). The same 
model presents also evidence of different degrees of clusterization at different heights and 
this is in contrast with the experimental observation 0. The model with gravity in both 
directions and without vibrating wall shows a stationary flow in the horizontal direction, 
where there are periodic boundary conditions: the profiles of the number density n{y), 
the X component of the velocity Vx{y) and the granular temperature T{y) as functions of 
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the distance from the bottom y are in very good agreement with the experimental profiles, 
showing a linear behavior in a broad region near the bottom which corresponds to the region 
where the collisions dominate the dynamics. This version of the model also shows strong 
evidences of density dependent clusterization and a non-Gaussian behavior near the bottom 
wall. The simplicity of the first setup has allowed us to exactly solve the hydrodynamics 
equations for n{y) and T{y) following the formulation of Brey et al. 0]: however it is not 



possible to obtain a matching condition between the bulk of the granular assembly and the 
vibrating wall which is responsible for the injection of energy. It seems to be an intrinsic 
problem of the high density gradients observed near the bottom wall which the Navier- 
Stokes approximation fails to describe. This suggests the need of a better description of the 
boundary layer, which should include higher order density and temperature gradients and 
also, at the level of the kinetics, the non-Gaussian velocity statistics and the effect of spatial 
correlations (clustering). 
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APPENDIX A: THE BIRD'S SCHEME FOR MONTE CARLO SOLUTION OF 

BOLTZMANN EQUATION 

The Bird's scheme, often called Direct Simulation Monte Carlo (DSMC), was designed 



in the 1960s ||T^ and its derivation was a priori independent of the Boltzmann equation. 
Recently its convergence to solutions of the Boltzmann equation in a suitable limit has been 
proved ||23|, reinterpreting it as a measure- valued stochastic process. 



The Bird's scheme can be formulated as a fixed time (At) step "molecular dynamics- 
like" simulation. At each time step the dynamics is separated in two distinct processes: the 
independent evolution of every particle and the collisions of near particles. The following 
algorithm (for the single time step) is the one we implemented in this work, which is a 
modification of the original scheme. 

• Free flow: each particle evolves independently following the equations of motion x= v, 
v= g with first order discretization. 

• Collisions: every particle i, during this time step has a probability pc of making a 
collision, so that Pc/At oc a (in fact the colhsion cross section is proportional to the 
diameter of the particles). If the particle collides, then another particle Xj, Vj is chosen 
with |xj — Xjl < tb (we call tb the "Bird radius", but it could also be thought as the 
"Boltzmann radius") with probability pij oc |vr|; for the pair i,j the postcoUisional 
velocities are calculated as they were at contact with a random choice of the collision 
parameter fi; this step is repeated for every particle. 

It is important to stress the fact that this is above all a Monte Carlo method to solve 
the Boltzmann Equation (^: in this sense microscopic (short range) details are lost, the N 
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particles themselves do not represent A^ real grains of the granular assembly but carry the 
space-time average information of many more particles. 



APPENDIX B: THE NUMERICAL COEFFICIENTS IN THE HYDRODYNAMIC 

EQUATIONS 

In section V the hydrodynamics of the first model is studied. The equations with the 



transport coefficients calculated by Brey et al. [^ are used. The coefficients needed in our 
case are the two thermal conductivities k and /i appearing in the expression of the heat flux 



q = —K\/{kBT) — fiVn 
and the coefficient ( of the dissipative term 



(43) 



-CksT. 



(44) 



In reference ||21[ the coefficients are given for the case d = 3 {d is the dimension of the 
space). We have taken the coefficients for d = 2 from an unpublished (to our knowledge) 



work of Brey et al. ||26[ and we have put them in the following form: 



where 



and 



— K = A{r 
-/i = 5(r, 
-C = C{r) 



iksT) 



1/2 






an 



iksT) 



1/2 



m 



1/2 



A{r) = -Ki{r)Ko 
B{r) = -/ii(r)fi:o 
C{r) = -Ciir)/vo 



(45) 
(46) 
(47) 



(48) 
(49) 
(50) 
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no=^ (51) 



vr 

/^i = 2 X ,^,y^^ (53) 

l + ci(r) 
^^ ^ ..(r)-40(r) (^^) 

Ci = ^(l-r^)(l + ^ci(r)) (55) 

/19 15 1 \ 

i/i(r) = (1 + r) (^- - -r + ^^(14 - 6r)ci(r) J (56) 

^n^)-32^7_25^^30^2(i_^)- (57) 

The coefficients A{r), B{r), C{r) are plotted in the figure (|19|) . In the same figure are 
also presented the coefficients a{r) = {A{r) — B{r))/C{r) and /3(r) = {A{r) — ii?(r))/C(r) 
appearing in the equation (0) and, finally, the coefficients a'{r) = (a(r) + /3(r))/(2«(r)) 
and P'{r) = (l/(2«(r)))^/^ appearing in the solution (PT|). 
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FIGURES 




FIG. 1. A sketch of the first model where the granular assembly is driven by gravity plus a 
(periodically or stochastic) vibrating wall 




FIG. 2. A sketch of the second model where the only energy source is gravity, with components 
in both directions 
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FIG. 3. Snapshots of the model a) with stochastic wall at temperature Ty, = 50 and Ty^ = 250. 
The leftmost inset displays the time-averaged number density profile for both case. Values of other 
parameters: A^ = 500, N,„ f» 56, r = 0.7, r,„ = 0.7,(7^ = —1 
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FIG. 4. Distribution of rescaled horizontal velocities vj \fT^ for the model a) with stochastic 
wall at different temperatures T^ = 50, T^j = 100, T^ = 250. The other parameters are N = 5000, 
iV,„ « 180, r = 0.7, r,„ = n 7 n, = -i 
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FIG. 5. Distribution of horizontal velocities, for the model a) with stochastic wall, measured 

on stripes at different heights and rescaled by the average temperature at that height. The inset 

shows the normalized number density profile with the position of the chosen stripes. N = 5000, 

Ny, ^ 180, r = 0.7, r^ = 0.7, g^ = -1, T^ = 100 
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FIG. 6. Granular (dimensionless) temperature T/{g(,rB) versus dimensionless height y/r^ 
(above) and versus number density n (bottom) for the for the model a) with stochastic wall, 
with A^ = 5000, N,fu ~ 180, r = 0.7, r^ = 0.7, g^ = —1. The solid line is a power-law fit for T{n). 
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FIG. 7. Snapshot of the model a) with periodically vibrating wall (right) and time-averaged 
density profile (left) for the following choice of parameters: N = 500, A^^^, ~ 56, r = 0.5, r^„ = 0.7, 
5e = -1, U = 4007r, A.^ = 0.1 
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FIG. 8. Distributions of horizontal velocities for the model a) with periodically vibrating wall 
for two different values of inclination, that is 5fe = ~1 and ge = —100, while the other parameters 
are fixed: N = 500, N^ « 56, r = 0.5, r^ = 0.7, U = 4007r, A^ = 0.1 
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FIG. 9. Cumulated correlation function C(R), as defined in the text, measured along stripes 
at different heights for the model a), with periodically vibrating wall. In the inset is displayed the 
number density profile, with the position of the chosen stripes. Here N = 500, Nyj ~ 56, r = 0.5, 



0.7, U = 4007r, A^ = 0.1 and g^ 
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FIG. 10. Normalized number density n, dimensionless horizontal velocity v^/ y/Ox'^'B a-^d dimen- 
sionless granular temperature T/ ^g^rs versus dimensionless height y/rs for the two dimensional 



inclined channel (model b): A^ 

(f, = 7r/6), r = 0.95, r^ = 0.95 
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FIG. 11. Normalized number density n, dimensionless horizontal velocity Vx/^Qx^B and dimen- 
sionless granular temperature T/^g^rs versus dimensionless height y/rs for the two dimensional 
inclined channel (model b): A^ = 500, A'^ ~ 56, gx = 1, gy = —2 (i.e.: the inclination angle 
(f, = 7r/6), r = 0.95, r^ = 0.4 
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FIG. 12. Cooling rate, as defined in the text, versus dimensionless height y/vB for the two 
dimensional inclined channel (model b): N = 500, A^^ ~ 56, Qx = 1, gy = —2 (i.e.: the inclination 
angle (p = ^/6), r = 0.95, r^ = 0.95 or r^ = 0.4 
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FIG. 13. Distribution of horizontal velocities for the model b), measured on stripes at different 
heights and rescaled in order to have the same mean and variance. The inset shows the normalized 
number density profile with the position of the chosen stripes. N = 500, N^ ~ 56, r = 0.95, 
r^i = 0.95, Qx = 1, gy = —2 (i.e.: the inclination angle (j) = tt/6) 
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FIG. 14. Cumulated correlation function C(R), as defined in the text, measured along stripes 
at different heights for the model b). In the inset is displayed the normalized number density profile 
with the position of the choose stripes. Here N = 500, N^ « 56, r = 0.95, r^ = 0.95, g^ = 1, 
Qy = —2 (i.e.: the inclination angle (j) = tt/G). The dashed lines represent the power-law fits, the 
vertical dot-dashed line represent the width of the stripes Ay 
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FIG. 15. Normalized number density n, dimensionless horizontal velocit y Vx/y/gxfB and di- 
mensionless granular temperature T/^g^rB) versus dimensionless height y/rs for the two dimen- 
sional inclined channel. Here tangential restitution coefficients smaller than one are considered (see 



text): N = 500, Nyj « 56, gx = 1, gy = —2 (i.e.: the inclination angle (p = ^/6), r" = 0.95, 
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FIG. 16. Plot of H, defined in section V, versus I, for three different simulations of the model 
a): two cases are with the stochastic wall {N = 5000, N^u w 180, r = 0.7, r^ = 0.7, ge = —1, 
Tuj = 150 and T^ = 250), while the third case is with the periodic wall {N = 5000, N^ ~ 180, 
r = 0.7, r^ = 0.7, ge = -1, U = SOvr, A^ = 0.1 
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FIG. 17. Profiles of dimensionless hydrodynamic fields n, Vy and T versus the dimensionless 
height y/rB^ for the model a) with the stochastic wall at temperature T^ = 250. N = 5000, 
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180, r = 0.7, r^ = 0.7, ge 



-1. The dashed vertical line marks the same height of figure 
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FIG. 18. Profiles of dimensionless liydrodynamic fields n and T versus 1 — l/lmax (the new 
coordinate / is defined in section V and Imax = <^/Lx ~ vb/Lx), for the model a) with the 
stochastic wall at temperature T^ = 250. N = 5000, N^ ~ 180, r = 0.7, r^ = 0.7, Qe = —1- The 
solid lines are the theoretical fit using the hydrodynamics model of Brey et al. The vertical dashed 



line marks the height (also appearing in figure (17)) where T presents a minimum and, therefore, 
goes to 0. 










9 




0.4 


0.6 


0.8 1 


1 


- 










U) 




— 


p 


10 




1 




^--^ 


1 


, 



0.2 0.4 0.6 0.8 




FIG. 19. Transport coefficients A, B and dissipative coefficient C of hydrodynamics (section 
V), numerical coefficients a and j3 of the equation (|30|), numerical coefficients a' and /?' of the 



solution pi] ) versus the restitution coefficient r. 
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